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Distributed Relaxation for Conservative 

Discretizations 
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A multigrid method is defined as having textbook multigrid efficiency (TME) if the so- 
lutions to the governing system of equations are attained in a computational work that is a 
small (less than 10) multiple of the operation count in one target-grid residual evaluation. 
The way to achieve this efficiency is the distributed relaxation approach. TME solvers 
employing distributed relaxation have already been demonstrated for nonconservative For- 
mulations of high-Reynolds-number viscous incompressible and subsonic compressible 
flow regimes. The purpose of this paper is to provide foundations for applications of 
distributed relaxation to conservative discretizations. A direc t c orrespon den ce be tween 
the primitive variable interpolations for calculating fluxes in conservative finite-volume 
discretizations and stencils of the discretized derivatives in the nonconservative formula- 
tion has been established. Based on this correspondence, one can arrive at a conservative 
discretization which is very efficiently solved with a nonconservatiye relaxation scheme 
and this is demonstrated for conservative discretization of the quasi one-dimensional Eu- 
ler equations. Formulations for both staggered and collocated grid arrangements are 
considered and extensions of the general procedure to multiple dimensions are discussed. 


Introduction 

Full multigrid (FMG) algorithms 1,2 ’ 9 ’ 18,22,23 are 
the fastest solvers for elliptic problems. These algo- 
rithms can solve a general discretized elliptic problem 
to the discretization accuracy in a computational work 
that is a small (less than 10) multiple of the opera- 
tion count in one target-grid residual evaluation. Such 
efficiency is known as textbook multigrid efficiency 
(TME). 3 The difficulties associated with extending 
TME for solution of the Reynolds-averaged Navier- 
Stokes (RANS) equations relate to the fact that the 
RANS equations are a system of coupled nonlinear 
equations that is not, even for subsonic Mach num- 
bers, fully elliptic, but contain hyperbolic partitions. 
TME for the RANS simulations can be achieved if 
the different factors contributing to the system could 
be separated and treated optimally, e.g., by multigrid 
for elliptic factors and by downstream marching for 
hyperbolic factors. The way to separate the factors 
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is the distributed relaxation approach proposed ear- 
lier. 1 ’ 2 Usually, distributed relaxation can be applied 
throughout the entire domain having the full effect 
away from discontinuities (shocks, slip lines) in the reg- 
ular (smoothly varying) flow field. Some local relax- 
ation sweeps should be applied in these special regions 
after (and perhaps before) the distributed relaxation 
pass to reduce residuals to the average level charac- 
terizing the regular flow field. The general framework 
for achieving TME in large-scale computational fluid 
dynamics (CFD) applications has been discussed. 6,21 

The approach to the solution of the RANS equa- 
tions proposed in this paper is based on an FMG 
algorithm with multigrid cycles employing distributed 
relaxation. It is envisioned that the FMG-1 algorithm 
(an algorithm with one multigrid cycle per level) will 
provide solutions with algebraic error below the level of 
the discretization error. Another useful characteristic 
of the solution process is a possibility to rapidly con- 
verge residuals to machine zero. The latter property 
is not necessary for achieving TME, but it is highly 
favored in practical applications. 

The distributed relaxation approach relies on a prin- 
cipal linearization of the governing system of nonlin- 
ear equations. The principal linearization of a scalar 
equation contains the terms that make a major contri- 
bution to the residual per a unit change. The principal 
terms thus generally depend on the scale, or mesh 
size, of interest. For example, the discretized high- 
est derivative terms are principal on grids with small 
enough mesh size. For a discretized system of differ- 
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ential equations, the principal terms are those that 
contribute to the principal terms of the determinant 
of the matrix operator. 

Development of a solution method for the govern- 
ing RANS system can be significantly simplified if 
the target discretization possesses two properties: ( 1 ) 
The discretization of the corresponding principally lin- 
earized system is factorizable , 2, 3,16,1 7 i.e., the determi- 
nant of the principal part of the discrete matrix oper- 
ator can be represented as a product of discrete scalar 
factors, each of them approximating a corresponding 
factor of the differential matrix operator determinant. 
(2) The obtained scalar factor discretizations should 
reflect the physical anisotropies and be efficiently solv- 
able. These two properties are not necessarily required 
for convergence of the algebraic errors in the FMG-1 
algorithm, but they are very important for ensuring a 
fast residual convergence and providing accurate solu- 
tions. 

Appropriate choices of factorizable discretizations 
for nonconservative formulations and textbook ef- 
ficient multigrid solvers employing distributed re- 
laxation have already been demonstrated for high- 
Reynolds-number viscous incompressible 7, 20 and sub- 
sonic compressible 19 flow regimes. An example of 
the distributed relaxation involved in calculation of 
the Euler flow with shock has been demonstrated ; 21 
the finite- volume collocated-grid discretization scheme 
used in 21 was a standard flux-differencing splitting 
scheme of Roe . 15 Applied to a one-dimensional prob- 
lem, the scheme is factorizable and provides reasonable 
approximations for the determinant factors. However, 
in multiple dimensions, the scheme is not factorizable, 
and other (factorizable) schemes should be considered. 
The purpose of this paper is to provide foundations 
for applying distributed relaxation to conservative dis- 
cretizations of the Euler equations corresponding to 
factorizable discrete schemes. 

The present material is organized in the following 
order: First, the Euler equations for inviscid compress- 
ible flow problems are defined. Secondly, the idea of 
distributed relaxation is briefly explained from both a 
differential and a discrete viewpoint. The attributes of 
a desired nonconservative Euler-system discretization 
scheme are discussed and a model problem, the one- 
dimensional Euler equations, is used to illustrate the 
derivation of conservative discretizations correspond- 
ing to a given nonconservative schemes. Numerical 
tests are reported for a collocated grid scheme for sub- 
sonic quasi-one-dimensional flow. Finally, concluding 
remarks are given. 

Euler Equations 

The time-dependent three-dimensional Euler system 
of compressible flow equations can be written as 

<9fQ-bR=0, (1) 


where the conserved variables are. Q = 
(pu, pv, pw, p, pE) T , representing the momentum 
vector, density, and total energy per unit volume, and 

R(Q) = d x F(Q) + 3 y G(Q) + cbH(Q), (2) 


F(Q) = 


pu 2 + p \ 


/ puv \ 

p uv 


pv 2 + p 

puw 

, G(Q) = 

pvw 

pu 


pv 

mE + up ) 


\ pvE T vp j 


H(Q) = 


/ puw \ 
pvw 

pw 2 + p 
pw 

\ pwE + wp J 


( 3 ) 


In general, the simplest form of the differential equa- 
tions corresponds to nonconservative equations ex- 
pressed in primitive variables, here taken as the set 
composed of velocity, pressure, and internal energy, 
q — (u 5 v, w, p, e) T . The primitive variables are related 
through the equation of state, 


(7 - !)/*> 

U) 

E - t (« 2 + v 2 + to 2 ) , 

(5) 

7 P/P, 

( 6 ) 


and 7 is the ratio of specific heats. 

The time-dependent nonconservative equations are 
found readily by transforming the time-dependent con- 
servative equations. 


fa[<9 (Q + R] = o, 
dtq + fgR = 0, 


where is the Jacobian matrix of the transforma- 
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For steady-state equations, the time derivative is 
dropped. 

In an iterative procedure, the correction <^q == 
qrc-H _ q n ? w h ere n i s an iteration counter, can be 
computed from the equation 


L *o = -f5 R ’ 


( 9 ) 


where L is the principal linearization of the noncon- 
servative operator. 


L = 


Q 0 0 fa 0 

0 Q 0 fa 0 

0 0 Q fa 0 

pc 2 d x pc 2 d y pc 2 d z Q 0 

C ^d x C ^8y 0 Q 


( 10 ) 


where Q — ud x T vd y + wd z = (u * V), and the coef- 
ficients u = (u, and c 2 are evaluated from the 

previous approximation q n and, for the current itera- 
tion, are considered as known values unrelated to the 
target primitive variables. The last equation in (9) is 
decoupled from the first four, representing the convec- 
tion of entropy along a streamline. The determinant 
of the matrix operator L is 


Q 3 [Q 2 -c 2 A], (11) 

where A is the Laplace operator, and Q 2 — c 2 A repre- 
sents the full-potential operator. Note that the right 
side of the correction equation (9) is a combination of 
conservative-discretization residuals. The left side ap- 
proximates the nonconservative equations. Thus, away 
from discontinuities in the flow field, we expect a good 
correction to the previous solution approximation. 


Distributed Relaxation: Differential 
Equations 

The distributed relaxation method for the Euler 
equations replaces <5q in (9) by M5w, 


M 


1 

0 

0 

0 

0 


Q 

0 


(12) 


so that the resulting matrix L M becomes lower trian- 
gular, as 


LM = 


Q 0 0 

o <5 o 

0 0 Q 

pc 2 d x pc 2 d y pc 2 d z 


0 0 ' 

0 0 

0 0 , 

Q 2 - c 2 A 0 
0 Q 

(13) 


LM<hv = -^R. (14) 

The diagonal elements of L M are composed of the 
factors of the matrix L determinant and represent the 
elliptic or hyperbolic partitions of the equation. The 
5 w variables were termed earlier, 2,8 “ghost variables,” 
because they need not explicitly appear in the calcu- 
lations. 

The distributed relaxation approach yields fast con- 
vergence for both steady and unsteady simulations, if 
the constituent scalar diagonal operators in L M are 
solved with fast methods. The approach can be ap- 
plied to quite general equations; a set of matrices M 
has been derived 2,3 ) that provide a convenient lower 
triangular form for the compressible and incompress- 
ible equations of fluid dynamics (including a variable 
equation of state). 

For the compressible Euler equations, the scalar 
factors constituting the main diagonal of L M are 
convection and full-potential operators. An efficient 
solver for the former can be based on downstream 
marching, with additional special procedures for re- 
circulating flows; 7,8, 14,24 the latter is a variable type 
operator, and its solution requires different procedures 
in subsonic, transonic, and supersonic regions. In deep 
subsonic regions, the full-potential operator is uni- 
formly elliptic and therefore standard multigrid meth- 
ods yield optimal efficiency. When the Mach number 
approaches unity, the operator becomes increasingly 
anisotropic and, because some smooth characteristic 
error components cannot be approximated adequately 
on coarse grids, classical multigrid methods severely 
degrade. The characteristic components are those 
components that are much smoother in the charac- 
teristic directions than in other directions. 3, 12,13 In 
the deep supersonic regions, the full-potential opera- 
tor is uniformly hyperbolic with the stream direction 
serving as the time-like direction. In this region, an 
efficient solver can be obtained with a downstream 
marching method. However, downstream marching 
becomes problematic for the Mach number dropping 
towards unity, because the Courant number associated 
with this method becomes large. Thus, a special pro- 
cedure is required to provide an efficient solution for 
transonic regions. A possible local procedure 4,5, 10 is 
based on piecewise semicoarsening and some rules for 
adding dissipation at the coarse-grid levels. 

Boundaries introduce some additional complexity in 
distributed relaxation. The determinant of L M is 
usually higher order than the determinant of L. In 
this case, it is higher by the factor Q. Thus, as a 
set of new variables, S w would generally need addi- 
tional boundary conditions. In relaxation, because the 
ghost variables can be added in the external part of the 
domain, it is usually possible to determine suitable 
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boundary conditions for <5w that satisfy the original 
boundary conditions for the primitive variables. Ex- 
amples are given in 20 for incompressible flow with en- 
tering and no-slip boundaries. However, to construct 
such a remedy may be difficult and/or time-consuming 
in general. In addition, enforcing these boundary con- 
ditions causes the relaxation equations to be coupled 
near the boundaries, not decoupled as they are in the 
interior of the domain. Thus, near boundaries and 
discontinuities, the general approach 2,3 is to relax the 
governing equations directly in terms of primitive vari- 
ables. Several sweeps of robust (but possibly slowly 
converging) relaxation, such as Newton-Kacmarcz re- 
laxation, can be made in this region. The additional 
sweeps will not affect the overall complexity because 
the number of boundary and/or discontinuity points 
is usually negligible in comparison with the number of 
interior points. 

Distributed Relaxation: Discrete 
Equations 

The discrete equations are formulated in terms of 
discretized primitive variables. One way to attain 
a conservative finite-volume discretization, R h , cor- 
responding to the discrete version of (1) requires an 
interpolation of the primitive variables to the inter- 
face locations, maintaining the required order of accu- 
racy. In this paper, the considered grid is a Cartesian 
structured grid, although generalization to curvilinear 
coordinate systems is possible. At this moment, we do 
not make any assumptions about locations of the vari- 
ables and volume interfaces, i.e., these locations can 
be staggered or collocated. 

The equations for the relaxation update are formu- 
lated similarly to (9) 

L h <Sq h = "lq Rh (15) 

The right side of the equations are a combination of 
the conservative residuals with the coefficients of Jq 
calculated at the points where the corresponding con- 
servative equations are defined. 

The derivation of the discrete operator L h requires 
special considerations. 

• The linear discrete operator L h should relate to 
R h as 

L” = < 16 > 

where 7 Z h is the principal linearization operator 
of R h . 

• L h should be factorizable. 

• To reflect the correct domain of influence, the dis- 
crete factors appearing at the main diagonal of the 


matrix product L h M h , where M h is the discrete 
distribution matrix (compare to (13)), should be 
type-dependent, i.e., central for elliptic factors 
and upwind (or upwind biased) for hyperbolic fac- 
tors. For Euler equations, it is also required that 
the discretization of the full-potential factor re- 
flects the physical anisotropies of the differential 
full-potential operator. 

Some simplification can be achieved by replacing L h 
with a less accurate operator providing a good conver- 
gence with defect-correction iterations. For staggered- 
grid nonconservative formulations, proper operators 
L h have already been derived for incompressible and 
low-Mach-number flows. 1 ' 2 ’ 19,20 For collocated-grid 
formulations, a family of factorizable discrete opera- 
tors L h is analyzed in. 11 The collocated-grid scheme 
presented in this paper belongs to this family. An 
alternative approach also resulting to factorizable dis- 
cretizations is proposed by Sidilkover in 16,17 

The basic scheme, L h basicj is similarly formulated 
for both staggered and collocated variable arrange- 
ments: 



' Q h 

0 

0 


o ' 



0 

Q h 

0 

P z 

0 


T h _ 

u basic — 

0 

0 

Q h 

0 

) 


pc 2 dx 

(>C 2 0y 

pc 2 d h z 

Q h 

0 



c±d h 

TT * 

c -d h 

7 v 


0 

Q h _ 

(17) 


where the discrete derivatives, , <9 * , in all off- 

diagonal terms are the second-order-accurate central- 
differencing approximation (short, i.e., with h - spacing 
for staggered grids and long, i.e., with 2/i-spacing, for 
collocated grids). All the diagonal terms, Q h , ex- 
cept Q h in the fourth equation, are discretized with 
the same second-order- accurate upwind (or upwind- 
biased) discretization scheme. In the subsonic regime 
(|u| 2 = u 2 + v 2 + w 2 < c 2 ), the (5^-term is differenced 
with a second-order-accurate downwind (or downwind- 
biased) discretization. 

The determinant of the matrix operator L h basic ^ 
given by 

(e A ) 3 [Q' , Q' , -c 2 A A ], (18) 

where A is a discrete Laplace operator. 

For staggered grids, the discrete full-potential op- 
erator appeared in the brackets is a reasonable h - 
elliptic operator. However, its stencil is somewhat 
wide (because of Q h Q h ) and does not reflect, the phys- 
ical anisotropies of the differential full-potential oper- 
ator. The collocated-grid full-potential operator also 
has a wide stencil but suffers a more serious drawback: 
The Laplace operator A* is a wide (with mesh spac- 
ing 2 h) operator. For slow velocities (|u| & 0), the 
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discrete full-potential operator is dominated by this 
non-/i-elliptic wide Laplacian; and efficiency of any lo- 
cal relaxation in error smoothing severely degrades. 

A way to improve the discrete full-potential operator 
is to change the discretization of Q h to Q h + A h . Then 
the discrete full-potential operator is changed to 

Q\4 a + Q h Q h -c 2 A h . (19) 

If the operator A h is at least second-order small (pro- 
portional to /i 2 ), the overall second-order discretiza- 
tion accuracy is not compromised. The ideal choice for 

A h is A h - (Q h y l V h ,V h = T h - {Q h Q h - c 2 A h ), 

where T h is a desired discrete approximation for the 
full-potential factor. We omit the discussion on what 
would be the best discretization for the multidimen- 
sional subsonic full-potential operator. Note only that 
it is possible to construct a discretization that satisfies 
the following properties: (1) At low Mach numbers, 
the discretization is dominated by the standard (with 
mesh spacing h ) ^-elliptic Laplacian. (2) For the Mach 
number approaching unity, the discretization tends 
to the optimal discretization for the sonic-flow full- 
potential operator (see 4,5,10 ). 

The operator is a nonlocal operator and 

its introduction can be effected through a new auxil- 
iary variable ip h and a new discrete equation Q h xj) h — 
V h p h . For the purpose of constructing a correspond- 
ing conservative scheme, the new variable should 
be further replaced with ip h — Q where Q* 
is an arbitrary discrete approximation to the oper- 
ator Q and <p h is another auxiliary variable. Thus, 
the new vector of discrete unknowns becomes q h = 


(iA v h 

,«AvA 


) T ■ 

The 

discrete 

operator 

L basic 

is changed to L h 







" Q h 
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0 
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0 
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0 
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1 

Q h 

0 

0 


0 

0 

0 

-1 

0 

<2* 

0 


cl d h 

L 

c -d h 

7 y 

-d h z 

7 z 

0 

0 

0 

Q h 


( 20 ) 

The corresponding distribution matrix, M h , for dis- 


tributed relaxation 

is defined 

as 
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0 
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-K 

0 
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0 

0 

1 
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0 

0 
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0 
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0 
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0 
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0 

0 
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0 

Q h 
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0 
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0 

0 
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so that the resulting matrix L h M h becomes lower tri- 
angular as 


Q h 

0 

0 

0 

0 

0 

0 

0 

Q h 

0 

0 

0 

0 

0 

0 

0 

Q h 

0 

0 

0 

0 

0 

0 

0 

Q h 

0 

0 

0 

pc 2 d h x 

pC 2 d h y 

pc 2 d z 

1 

T h 

0 

0 

0 

0 

0 

-1 

0 

Q* 

0 

c -d h 

7 U * 

cld h 

7 u v 

: c -d h 

7 2 

0 

0 

0 

Q h 


( 22 ) 

In the following sections, a direct correspondence 
between the primitive variable interpolations for cal- 
culating fluxes in the conservative discretization and 
stencils of the discretized derivatives in the noncon- 
servative formulation is established. This correspon- 
dence can be used in several ways. For example, one 
can define interpolations of variables in a conserva- 
tive discretization so that the corresponding operator 
L h from (16) coincides with a predefined linear oper- 
ator derived from a nonconservative formulation. As 
another example, one can derive the linear operator 
L h corresponding to a given conservative discretiza- 
tion and then try to find an appropriate matrix M h to 
design an efficient distributed relaxation. Some mixed 
requirements, partially defined by the interpolations 
in the conservative discretization and partially by the 
discretization stencils in L h , are also possible. 

A brief description of the derivation of the 
conservative-nonconservative correspondence can be 
done as follows. The starting point is the full Newton 
linearization of the conservative discretization. This 
linear operator acts on the perturbation function. The 
assumption of smoothness suggests that the changes 
of primitive variables on the scale of the characteris- 
tic meshsize are small in comparison with the absolute 
value of the function. Under this assumption the prin- 
cipal linearization, H h , retains only terms including 
the perturbation function derivatives rather than the 
perturbation function itself. The product l\ h ap- 
proximates the set of nonconservative equations; and 
one can use it to establish the relations between the 
flux interpolations in the conservative discretization 
and the discrete differencing in L h . 

In addition to regions with shocks and slip lines, 
the smoothness assumption is not valid in the regions 
where primitive variables have small absolute values, 
e.g., in the neighborhood of a stagnation point; in 
these regions, therefore, distributed relaxation should 
be replaced with local relaxation. 

One-Dimensional Euler System 

The steady-state one-dimensional equations express 
the conservation of momentum, mass, and total energy 
as 
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where 


R-(Q) = 0, 


Q = (pu, p, pE) T , 

R = d x (F). 


(23) 


(24) 

(25) 


The flux F is defined as 



(26) 


The vector of primitive variables is q = (u,p,e) T . 
The Jacobian matrix of the conservative-to- 
nonconservative transformation is defined as 


5q 

dQ 


The set of the nonconservative equations is defined 
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where the capital letters designate the interpolations of 
the primitive variables to the flux locations ( E repre- 
senting interpolation of e). Note that these interpola- 
tions can be different for different equations, therefore, 
breve, bar, and hat notations for momentum, mass, 
and energy conservation equations respectively. Under 
the smoothness assumption, the principal linearization 
7 Z h of R h applied to disturbance function, S q, results 
in a set of the three linear operators. 


a h (s q ) = 


2 PU (frQright ~ (‘Heft 
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\ Wright - (Heft 
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The discrete matrix operator corresponding to (28) is 
defined as 
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uD? 

Ulf)P 
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uD e 


\ 


/ 


(29) 


where the D-entries are some discrete operators ap- 
proximating the first derivatives from (28); the 're- 
entries approximate zeros, meaning that the difference 
operators can either be zeros or approximate a higher 
order derivative multiplied by a power of h] the su- 
perscript denotes the variable to which the difference 
operator is applied. In most cases, the operators are 
computed at integer locations except for staggered ar- 
rangements, where the breve operators are defined at 
half-integer locations. 

The conservative discretization of R (23) can be ex- 
plicitly written out as 
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where the lZ h operators are defined at the same lo- 
cations as the corresponding operators in (28), i.e., 
for collocated grids - at integer points, for staggered 
grids - the first operator, it h , is defined at half-integer 
points. Tilde values are arbitrary interpolations satis- 
fying the desired accuracy order, and the A-notation 
is self-explanatory. 

Using (27), an approximation to the first equation 
of (28) is obtained as 

( 7 zjll 1 (ii h - un h ^j = u (2A U - A u ) 

+ ilzi^ AP 

P 

+j ( A ' ‘ A ') 
-^(A‘-A‘) (31) 

For staggered grids, the bar-operators (7J ft ,A u ,A p , 
and A € ) are averaged from the neighboring locations. 
Approximations to the second and third equations of 
(28) are similarly computed: 

( 7 - 1 ) (~-un h + ^-P h + n h ^j = 
yPA u + ^ (— 2A U + \A U + §A U ) 

+U (yA p - (7 - 1)A P ) + ^ (— A p + |A p + |A p ) 

(a c - IV - 1 V) , 

(32) 

(- un h + (-£ + n h + n h ) = 

E (— A“ + 7A u ) + U 2 (-2A U + IA !1 + f A u ) 

+ {2^m (_A p - ^Ap + ^Ap) 

+^(-A p + iAp + iA p ) 

+UA C + ^ (A f - IA ( - |A e ) . 

(33) 


For staggered-grid formulation, breve- values are 
computed as averages of corresponding values at the 
half- integer points. 

Comparison of (31)— (33) with (29) gives the follow- 
ing relations between conservative and nonconserva- 
tive discretizations. The recurrent representation form 
is chosen for the sake of compactness. 


A u =~fD u -( 7 - 1)£> U , (34) 

A u = \D" + ±A“, (35) 

A " = (d u - /? ^— 2 A u + l - A u ) ) , (36) 

A P = D P -V P , (37) 

,38) 


\7 V 7 V "2 )' ' 

A e = D e - V e , (40) 

A c = -p-V< + A e , (41) 

7/? 

(42) 


where M 2 = u 2 /c 2 = u 2 /( 7(7 — l)c) is the local 
(squared) Mach number at the point defined by the 
subscript index j, fi = (7 — 1 ) A / 2 , 9 = y+Vj, an< ^ 


Derivation of the conservative discretization 
corresponding to a given nonconservative scheme 

Let us consider three different discretizations of (29) 
associated with the subsonic flow regime. In all the 
cases the P-entries are zeros and other off-diagonal 
terms are second-order accurate central discretizations 
of the first derivative. 

A First-order scheme: D u and D € are the first-order 
upwind discretizations, D p is the first-order down- 
wind discretization. 

B Second-order biased scheme: D u and DJ are the 
second-order upwind discretizations, D p is the 
second-order downwind discretization. 

C Second-order central scheme: D u , , and D ^ 
are the second-order central discretizations. 

The following formulas define the interpolations for 
evaluating fluxes in the conservative discretization 
(30); the superscript NC denotes the fluxes in rep- 
resentation of the D-derivative terms of (29) (see Ta- 
ble 1 .) 
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Table 1 The fluxes in representation of the /^-derivative terms of (29). 



A 

B 

c 

r\u ->+2 J 2 

f)U Vj — Vj-i 

JJU %Uj-2uj-i + ^Uj-2 


D i - — h — 

V J - h 

+ \ — \ U j - 

v ~ 2 h 

U "+ 1 = f“i+l + 

pNTT _pN'C 

2 3 ~ 2 
~ h 

f)P _ Pi + i-Pi-i 
^3 ~ 2 h 

P/4 - fa+\ + fa 

r\P — Pi+i-Pi-i 

U S ~ 2 h 

= fa +1 + fa 

nP — P>fi~Pi-=i 
v “ 2h 

Pj + l ~ \ Pj + 1 

n« 2 3 2 



f)U _ ' U 7 + 1“ U J-1 

D i - — h — 

U J ” 2 h 

= 2 u i + 1 + 2 U J 

_ u i ~ , 2 ft 

U N + i = fa+i + \uj 

~ 2 h 

Ufo = 5“i+i + 

pNO y'p.VT — 

rSP J+2 3 ~ i 

f)P _ Pl + l-Pj 

jjp _ - iPi+2 + 2 Pj + i“TPj 

f\P Pj + 1 —Pj~l 

D j - — h — 

~ h 

pNC — n . , 

Pfal = |pj + l - fa + 2 

U J ~ 2 h 

pjii = & + fa+i 

T Ml 2 3 2 

j\u — u i+x ~ u i-i 

nu u i+i~ u i-\ 

n« — u j+ i r:5 t j- 1 

°J ~ h 

U j “ 2 h 

Uf+i = fa + 1 + fa 

U j ~ 2h 

U? + C i = 

— 2h 

&j+i = 2 u j+ x + 

fSe •*+ 2 2 

r\e € i~ 


f)e _ f j+i 

D j ~ h 

u i ~~ h 

U i - h 

j?NC _ 3 , 1 . . , 

&j + l ~ 2 € J 

U 3 “ 2/x 

p.VC' _ 1 r . 1 , . 

^j + 1 “ 2 e J 2 G + 1 



= + 


(43) 


= 5^4 + 

5 ^' 


(44) 

^ + i 

= 1? 

- 

- ;) (-2C 

i+i + 2 ^+ 5 ) 

(45) 


__ p^c 



(46) 


1 

“ i+ 7 m 2 

(^4 + 

7A/ 2 F J+i ) , 

(47) 


= c 

pNC 

_J + i 

■7^+1 



+C(t 

-.)(> 

+ 7^ +i 

P 3 + k) 

(48) 


- E\ 

iVC 

i+e 



(49) 


- E j+h' 



(50) 


= 



(51) 


where C = 

Details of a collocated-grid scheme 

The flux computations for a collocated-grid scheme 
can be constructed more simply by separating the mass 
flux contributions in each equation. Using the defini- 
tions in Table 1, the flux contributions at the interface 
locations j + \ are given below. 


{P u ) 


pNCfjNC 

{-y-l)E NC ' 


(52) 

(53) 

[puE + pu) = ( pu)E NC + P NC U NC . (54) 


(pu 2 + p) = ( pu)U NC + P 


?NC 


The factor <j> is introduced into the mass flux so the 
scheme remains conservative, as 


<4 = 7/j - (55) 

where the fully-upwind scheme is used and Q* is taken 
as Q. Since <f) is second order small, zero conditions 
are imposed at the boundaries. 

Numerical Tests 

Computational results are shown in this section for 
the quasi-one-dimensional Euler equations. A uniform 
grid of N points is used over the domain 0 < x < 1 
with h — l/(N — 1). The area distribution is defined as 
<r(a:) = 1 — 0.8z(l - x). The flow is fully subsonic and 
the inflow Mach number is 0.5. For all of the results 
presented below, we overspecify values from the exact 
solution of the differential problem at the boundary 
and at locations outside the domain; hence conserva- 
tive fluxes are established at the half-grid locations 
surrounding the interior points. 

The maximum discretization errors versus grid size 
in Fig. 1 shows a second-order spatial convergence. 
The results were converged using a Full Multigrid 
(FMG) cycle starting from a converged solution for 
h = 1/16 and using two relaxations of a nonconser- 
vative operator, Eq. (15), at each subsequent level. 
In Table 2, the L 2 -norm of the discretization error 
in pressure, e d = pj - p(xj) exact i corresponds to the 
results with zero algebraic erors. The ratios of the al- 
gebraic error in this L 2 -norm to the discretization error 
in this norm for the FMG cycle are quite small at each 
grid. Because the nonconservative operator is factoriz- 
able and distributed relaxation has been demonstrated 
previously to be efficient solvers for nonconservative 
discretizations, we infer that distributed relaxation 
would yield high efficiency for this conservative dis- 
cretization. It remains to demonstrate this for more 
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h 


Fig. 1 Maximum discretization errors versus grid 
size for model problem. 


h 

IMI : P 

Ikall/lledll : p 

1/32 

0.7323xl0- 3 

0.043 

1/64 

0.1618xl0 _a_ 

0.086 

1/128 

0.3856xl0 -4 

0.006 

1/256 

0.9447xl0~ b 

0.004 


Table 2 The L 2 -norm of errors in p at convergence 
and the ratio of algebraic-to-discretization errors 
after an FMG cycle with two nonconservative re- 
laxations on each grid. 

general situations, including applications with more 
general boundary conditions and to flows with cap- 
tured shocks. 

Concluding Remarks 

This paper has provided some foundations for the 
application of distributed relaxation to conservative 
discretizations. A direct correspondence between the 
primitive variable interpolations for calculating fluxes 
in conservative finite-volume discretizations and sten- 
cils of the discretized derivatives in the nonconser- 
vative formulation has been established. Based on 
this correspondence, one can design a conservative 
discretization which is very efficiently solved with a 
nonconservative relaxation scheme. This has been 
demonstrated for a conservative discretization of the 
quasi one-dimensional Euler equations on a collocated 
grid. Formulations for both staggered and collocated 
grid arrangements as well as extensions of the general 
procedure to multiple dimensions have been discussed. 
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